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Abstract 

Recent studies of cluster distribution in various ecosystems revealed Pareto statistics for the 
size of spatial colonies. These results were supported by cellular automata simulations that yield 
robust criticality for endogenous pattern formation based on positive feedback. We show that this 
self-organized criticality is a manifestation of the law of proportionate effect, first discovered in the 
context of business firm size. Mapping the stochastic model to a Markov birth-death process, the 
transition rates are shown to scale linearly with cluster size. This mapping provides a connection 
between patch statistics and the dynamics of the ecosystem; the "first passage time" for different 
colonies emerges as a powerful tool that discriminates between endogenous and exogenous clustering 
mechanisms. Imminent catastrophic shifts (like desertification) manifest themselves in a drastic 
change of the stability properties of spatial colonies, as the chance of a cluster to disappear depends 
logarithmically, rather than linearly, on its size. 
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In spatially extended ecosystems the observed patchiness supplies a lot of information on a 
given system's parameters and may signal that the system is approaching a transition jl| . Of 
particular interest are ordered {2} 3, % 5] and disordered j(| vegetation patterns in the semi- 
arid climatic zone, where climate changes and overgrazing lead to an accelerated process 
of desertification [7|. Recently, two groups of researchers simultaneously published new 
results about patch-size statistics. Using high-resolution remote sensing techniques, these 
groups found power-law distributions of vegetation patch size (the probability of finding a 
vegetation patch of size n falls like where 3 is greater than unity and can vary between 
different regions) along, the Kalahary Transect , and in three Mediterranean ecosystems . 
Another new study 10| has shown that the spatial distribution of ant colonies in a tropical 
agroecosystem is also Pareto-like. Power-law statistics indicate the absence of an intrinsic 
scale in the system, and in many cases its appearance signals that a stochastic system is 
"at criticality" , i.e., experiencing a continuous phase transition ll|. It is surprising to find 
scale-free distributions along a wide range of environmental conditions; in most systems it 
appears only at the transition point [3]. The challenging problem, thus, is to identify the 
mechanism that yields robust critical behavior for different ecosystems under a wide scale 
of environmental conditions. 

To support and explain their findings, all three groups suggested (slightly different) cel- 
lular automata models that support self organized criticality. Effectively, there are only two 
parameters in these models: the global constraint on the overall fractional tree cover, /*, and 
the local strength of positive-feedback between plants (for the sake of concreteness we use 
hereon the spatial vegetation terminology; our results, however, are generic). Facilitation 
implies that the chance of an occupied cell to become empty decreases, and the chance of 
an empty cell to become occupied increases, if its neighboring cells are occupied. It turns 
out that these two constituents are sufficient to ensure self-organization and power-law dis- 
tribution of patch sizes, while the details of the model (like the definition of neighborhood, 
kernel function, etc.) affect only the slope 3. These results increase the need to identify the 
underlying feature that makes criticality a robust property of models with global capacity 
constraints and positive-feedback. 

In order to elucidate the dynamical law that governs self-organized criticality, we have 
mapped the cellular automata used by Scanlon et. al. into a Markov chain. The cellular 



automaton introduced by [8J was used here with only slight modifications. Every cell is either 
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empty (0) or covered by tree canopy (1), and local facilitation determines the transition 
probabilities between these two states. The actual overall fraction of occupied cells, /, 
affects the dynamics if it deviates from the fractional tree cover allowed by the environmental 
conditions (the carrying capacity) /*. The transition rates in the cell are, 
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py(i-o) = PiJ +L-L, (i) 

where pij reveals the effect of facilitation by other biomass units in the local neighborhood, 

1 I 1 if site (k, I) is populated 

Pi,j = -q a k,h °k,i = \ (2) 
15 k,i£S I if site (k, I) is empty 

where S is the area of a region of diameter r centered at the cell. Grazing was added to 
the model by letting a fraction G of the "death" events occur at a constant rate, independent 
of the local density p it j. 

During the simulation, clusters were tracked and their sizes ware monitored. The birth 
and death rates for size n clusters were determined by dividing the number of transitions 
n — > n + 1 and n — > n — 1, respectively, by the average number of size n patches. The 
parameters 7, a and Ao were extracted from b n and d n using linear regression. 

Switching to patch size oriented dynamics, the events where a patch of size n is grew or 
decreased were monitored to yield the corresponding transition rates. It turns out that a 
stochastic birth-death process (where only ±1 transitions are allowed) is a remarkably good 
approximation for the real dynamics. If b n is the birth rate (the chance of a patch of n cells 
to grow by one cell, or one vegetation unit, to n + 1) and d n is the corresponding death rate, 
the time evolution of the patch size distribution P n (t) is described by the master equation: 

J t Pn{t) = -{b n + d n )P n {t) + d n+1 P n+1 (t) + &^iP n _i(t), (3) 

and its steady state ir n should satisfy the local balance equation ir n b n = ir n+ id n +i. This 



allows for a recursive solution for the steady state 13J, once the birth-death transition 
probabilities are given. 

Figure 1 shows the dependance of b n and d n on n. Apparently, both b n and d n grow 
linearly with n. However, the death rate is larger by a constant A, i.e., 

b n = in d n = in + A. (4) 
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FIG. 1: Tracing the stochastic simulation along time, the rates in which clusters of size n grow by 
one biomass unit, b n , or decay by one unit, d n , were monitored (see methods, here the parameters 
are /* = 0.1, r = 3). Neglecting merging and splitting of colonies, where size is changed by more 
than one unit, one gets the Markov birth-death process illustrated in the upper panel. In the lower 
panel the birth-death rates (in arbitrary units) are shown in the absence of grazing pressure; the 
rate of both processes grows linearly with the colony's size, where the death rate is always larger 
by a constant A. In the inset, the difference A(n) = d n — b n is presented without (a) and with 
(b) grazing pressure together with a linear fit. While A(n) is approximately constant without 
grazing, it shows, on average, linear growth where grazing is present. This gives the ground for 
the approximation A(n) = an + Ao, where a increases with grazing pressure. 

Using the local balance equation the steady state may be obtained analytically: 

where ttq is determined by the normalization condition. The product (jSJ) may be calculated 
exactly, and then approximated for n > 1 + — to yield, 



7 r(n + l + 
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where A is a normalization constant. Comparison between the distributions obtained from 
the cellular automata simulations and the power-law ([6]) using the calculated values of 7 
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and A shows excellent agreement (Figure 2, upper panel). This supports the approximation 
made by replacing the general Markov chain with a birth-death process. 

The monotonic increase of event rates with cluster size emerges, thus, as the dynamical 

n 

law that generates criticality. This, in fact, is the law of proportionate effect [14|, introduced 



many years ago in order to explain the Pareto distributions in economic systems [15l . [16, 



17j. Similar mechanisms were shown to underly surname statistics [181 ]. relative species 



abundance [20f] , urban population size [19(] and many other observations. Although here the 
geometric random walk is biased towards extinction, the multiplicative nature of the process 
guarantees power-law distributions, where the slope is determined by the bias A. The results 
presented here are for facilitation radius r = 3; testing the model under different conditions, 
the slope (3 was found to increase with r. Approximating the real Markov chain as a birth- 
death process turns out to be exceptionally good; for all cases tested by us (r = 2, 3, 4) the 
slope obtained from A and 7 differs only by 5 — 10% from the slope measured directly on 
the cellular automata. 

As indicated by the field data analyzed by Kefi et. al., patch distribution deviates from 
power-law in the presence of grazing pressure. A truncated power-law has been suggested in 
that case, where the strength of the grazing determines the typical patch size above which the 
distribution crosses over to exponential decay. Introducing grazing into our model through 
a density independent mortality rate (see methods), the same truncated power-law has been 
obtained. Mapping to a birth-death Markov chain, the only qualitative difference is that, 
under grazing, the leftward bias A grows linearly with the patch size, i.e., A = A + an (See 
inset of Figure 1). Substituting in ([5]) one attains, 



b ^ i»r(i + ^) / 7 

ir r , - 



( 7 + a) 2 r(l + n + ^) V« + 7/ 
In the limit of n 3> 1 this expression becomes the truncated power-law, 



(7) 
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An~y 1+ ^)e- n/no , (8) 
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where tiq = is the patch size above which n n decreases faster than a power law. 



The equilibrium distributions (Ej) and ([8]) fit perfectly the power law and truncated power 
law obtained from direct numerics, as seen in fig. 2. 

The main characteristic that distinguishes endogenous from exogenous cluster formation 
mechanisms is the dynamic of the system. If spatial heterogeneity plays a major role and 
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FIG. 2: The equilibrium probability n n is presented vs. n on a log-log scale, showing power-law 
behavior in the absence of grazing [upper panel, (a)]. The circles are results attained from the 
stochastic cellular automata, where the full line corresponds to the analytic expression for the 
solution of the birth-death process. The lower panel (b) shows the same in the presence of grazing 
{G = 0.6, see methods), where the truncated power law attained from the simulations is compared 
with . The relevant parameters a 7 and Ao are those extracted by linear regression from Figure 
1. 



determines the location and the size of large colonies, the chance of such a colony to disappear 



due to demographic stochasticity is exponentially small 



2l| , so in practice, on experimental 



;ime scales, large colonies never disappear. This feature was used by Vandermeer et. al. 
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who utilized the appearance/disappearance of large ant colonies as a manifestation of 



endogenous population dynamics. Our mapping of the system into a Markov chain allows 
us to quantify these characteristics and to establish a connection between the steady-state 
distribution and the dynamics of a single cluster. 



We define t™ as the mean first passage time 22j from some initial cluster size n to some 
smaller size m. Since the mean time for a cluster of size n to undergo a birth/death step is 
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FIG. 3: r n _ m , The average time elapsed until a colony of size n reaches a smaller size m at the 
first time is plotted against n — m for different environmental conditions. In the absence of grazing 
(G = 0, red circles) the stability of a colony grows linearly with its size. Under grazing pressure 
(G = 0.4, green and G = 0.6, blue) the time to extinction grows only logarithmically with the 
patch size. Datasets were taken directly from the simulation of the cellular automata with r = 3 
and /* = 0.1, and are fitted (full lines) to the predictions (110p using nonlinear curve fitting. 

l/(frn + d n ), the mean first passage times should satisfy: 



m _ 1 , b n d n 

' b n + dn b n + d n n+l b n + d n n 



Substituting b n = •yn and d n = jn + a n + A , we get, for n>l: 

^sK 5 ^)- (io) 

In that limit, the first passage time depends only on the difference n — m. In the absence of 
grazing (a — > 0), we get r n _ m = (n - m)/A . 

Evidentally, under grazing pressure, the dependence of r„_ m on n—m is logarithmic, while 
in the absence of grazing, the dependence is linear, as demonstrated in Figure 3. Tracking the 
temporal variations of colonies of different sizes (e.g., using satellite images taken yearly or 
seasonally) is thus an effective tool for detecting desertification. If large patches remain large, 
and their chance to reach a fixed, small size within some time frame decays exponentially 
with their size, it is plausible that the spatial structure reflects spatial heterogeneity, effects of 
local topography and so on. If the rate of disappearance for large clusters depends linearly 
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on size, it indicates that positive feedback plays a major role in the clustering process. 
Grazing, and other processes that lead to desertification, are characterized by logarithmic 
dependence on size differences. These three characteristic dynamics - exponential, linear, 
and logarithmic - differ strongly from each other, and thus temporal tracking provides a 
very reliable indication of the state of the system. 

It is interesting to note that, as firm size is governed by the same dynamical law, similar 
techniques are used for financial risk assessment. A Firm defaults when the value of its 
liabilities exceeds its assets. Merton [23| estimated the distance to default (DD) by modeling 



a firm assets as a geometric random walk, where the "step size" is proportional to the 
volatility. The expected default frequency 24| (EDF), which is the critical parameter that 
determines credit risk, is in fact an estimation of the typical time window needed for the 
geometric random walker to reach the default point; this quantity is just the inverse of the 
first passage time considered here. 



The models of 



flfl 



lOj all use a fixed fraction of occupied cells /*, i.e., a fixed total 



population. The differences in colony statistics quantify, in fact, the relative importance 
of the density-dependent growth factors (local facilitation, Allee effect, schooling) versus 
density-independent resources, like annual precipitation or the availability of food. If the 
patch distribution is Pareto-like with a relatively small slope j3, the system admits larger 
patches; this implies that the effect of positive feedback is strong enough to oppose stochas- 
ticity (purely stochastic dynamics lead to Poisson-like statistics). A truncated power-law, on 
the other hand, shows that density-independent factors control the overall carrying capacity. 
Given two systems with the same tree coverage /* but with different cluster statistics, the 
scale- free one is more robust against environmental changes (e.g., decreased annual rainfall) 
than the system characterized by a truncated power law. On the other hand, a strong Allee 
effect may result in irreversible transitions if the density falls below some critical level. The 
dynamical law presented here allows one to identify the transition between these states when 
the total population is still fixed; qualitative changes in the stability of large patches - from 
linear to logarithmic - may thus serve as an important precursor of catastrophic shifts 25] 
like desertification in arid ecosystems. 
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